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The problem of hydraulic fracture propagation is considered by using its recently suggested mod- 



ified formulation in terms of the particle velocity, the opening in the proper degree, appropriate 
spatial coordinates and e-regularization. We show that the formulation may serve for significant 
increasing the efficiency of numerical tracing the fracture propagation. Its advantages are illustrated 
by re-visiting the Nordgren problem. It is shown that the modified formulation facilitates (i) pos- 
sibility to have various stiffness of differential equations resulting after spatial discretization, (ii) 
obtaining highly accurate and stable numerical results with moderate computational effort, and (iii) 
sensitivity analysis. The exposition is extensively illustrated by numerical examples. 
G\ ' 

^ ■ 1 Introduction 

m 

Hydraulic fracturing is a widely used method serving to increase the linear size of an area of fluid or 
gas flow (see, e.g. the reviews in papers [TJ], [7J, [1]). In view of practical significance of the method, 
numerous papers have been published on the theory and numerical modeling of hydraulic fractures 
starting from the first publications [T8] , [I] , [22] , [II] , [H] , [21] , [2H] and [21] • The general formulation 
k> ' of the problem is well established (e.g. [2]). It includes (i) the fluid equations for flow of incompressible 

viscous fluid in the narrow channel; (ii) the solid mechanics (commonly static linear elasticity) equations 
defining the dependence of the channel's height on the pressure acting on the walls of the channel; (iii) 
equations of fracture mechanics defining the possibility of the fracture propagation and the trajectory of 
the fracture contour. Additional equations for proppant movement are added when accounting for the 
proppant injected at some stage of fracturing. 

The formulated mathematical problem is difficult from the computational point of view because 
of three major complicating factors: strong nonlinearity even in the simplest case of Newtonian fluid, 
caused by the fact that the channel's height (fracture opening) raised to some power enters the Poiscuillc 
equation as a multiplier by the unknown flux; moving boundaries of the fluid front and fracture contour; 
and, in general, the need to check the fracture conditions for finding the value and direction the frac- 
ture increment at time steps at each point of the front. Consequently, many investigations have been 



performed tending to reveal those general properties of the solution, which may serve for the problem 
simplification. They have provided knowledge on the asymptotics of the solution, possibility to neglect 
the lag between the fluid front and the fracture contour and on the typical regimes (e.g. [25] , Ej, Eh 
PP, [27], [6J, [3J, [5J, [23J, H3' [13, [IH])- The knowledge was incorporated in the computational codes 
for practical applications (e.g. [TBJ, [2]). Still, the mentioned difficulties are not overcome, and as 
emphasized in the review [2], there is the need "to dramatically speed up" simulators. 

The recent studies tended to address this challenge [19] - [22] have disclosed important features of 
the hydraulic fracture problem, which may be employed for enhancing the numerical simulation. They 
have led to the modified formulation of the problem using: 

1. the particle velocity, as a variable with continuous spatial derivative near the fluid front, instead 
of the pressure; 

2. the opening taken in a degree, defined by its asymptotic behavior at the fluid front, instead of the 
opening itself; 

3. the speed equation (SE) at each point of the front to trace the fracture propagation by the well- 
developed methods (see, e.g. [25]), instead of the commonly employed single equation of the global 
mass balance; the speed equation also presents the basis for proper rcgularization; 

4. e-rcgularization, that is imposing the boundary condition and the speed equation at a small distance 
from the front rather than on the front itself, to exclude deterioration of the solution near the front 
caused by the disclosed fact ([19], [20]) that the boundary value problem is ill-posed when neglecting 
the lag; 

5. the spatial coordinates moving with the front and evaluation of the temporal derivative under fixed 
values of these coordinates; 

6. reformulation of the common system of equations and boundary conditions in terms of the sug- 
gested variables complimented with e-regularization. 

The computational advantages of the modified formulation have been demonstrated [ID], [2D] by 
revisiting the classical Nordgren |25| problem. For it, common (without e-regularization) time stepping 
procedures could not provide reliable third digit and they led to strong deterioration of the solution 
near the fluid front for fine meshes. In contrast, applying e-regularization easily provided the solution 
with relative error less than 10~ 4 in a time stepping procedure; the solution was extremely stable and it 
never deteriorated near the front. Using of the suggested variables made it also possible [22] to obtain 
analytical solutions of the Nordgren [25] and Spence and Sharp [25] problems. 

In this paper, we make a further step in employing the modified formulation for (i) studying the 
stiffness of the system of differential equations arising after spatial discretization, (ii) increasing the 
efficiency of numerical tracing of the fracture propagation, and (iii) studying the sensitivity. To simplify 
the exposition and to compare the numerical results with benchmarks, we address the Nordgren problem. 

The structure of the paper is as follows. In Section 2, we briefly review various formulations of the 
problem. They include the common formulation (Subsection 2.1), the mentioned modified formulation 
(Subsection 2.2), its specification for the Nordgren problem (Subsection 2.3) and the self-similar formu- 
lation for ID problems. The latter serves us to write down the benchmark analytical solution of the 
Nordgren problem in the case of zero leak-off [22] and to obtain a benchmark solution for non-zero leak-off 
(Subsection 2.4). Further analysis employs the modified formulation and the benchmark solutions. 

Section 3 presents alternative approaches to spatial discretization of the lubrication partial differen- 
tial equation (PDE) and the speed equation. It is shown that they result in various systems of ordinary 
differential equations (ODE) with quite different stiffness. It appears that in schemes avoiding approxi- 
mation of the second spatial derivative by assuming the particle velocity fixed at an iteration step, the 



stiffness, in general, is not high being of order O(N), where N is the number of nodal points. In the 
case of constant particle velocity, the stiffness becomes of order 1/e independently on N (e is the reg- 
ularization parameter). This indicates favorable features of iterative schemes with the velocity fixed at 
the stage of integration in time. In contrast, in schemes employing approximation of the second spatial 
derivative, the stiffness depends on the number N more strongly: it is of order TV 2 for the Nordgren 
problem and of order iV 3 in the general case when the net pressure is connected with the opening by 
exact equations of the elasticity theory. 

Section 4 presents two alternative approaches based on the approximation of the second spatial 
derivative. The first of them employs the possibility to add the SE to the system of ODE, resulting from 
the spatial discretization of the modified lubrication equation. In this way, we obtain a new well-posed 
formulation for the joined system of ODE with initial (Cauchy) conditions on the fluid front. The 
formulation opens the possibility to employ methods, like those of Runge-Kutta, for solving the ODE. 
It serves us to use a standard MATLAB solver in further evaluations of the accuracy and sensitivity. 
The second approach employs classical Crank-Nicolson scheme to reduce the problem to tri-diagonal 
algebraic system. In this case, non-linear factors in front of the derivatives are iterated within a time 
step. The SE is used at iterations to find the new location of the fluid front. The rest of Section 
4 contains numerical results for these approaches and their discussion. It appears that the numerical 
procedures resulting from each of them are highly accurate, stable and robust. For comparison, we also 
present results obtained by using the equation of the global mass balance instead of the SE. It appears 
that the accuracy becomes notably (an order, at least) less than that when using the SE. 

In Section 5 we study sensitivity of the solution to changes of influx which is one of the major 
parameters of the problem. Conclusions are drawn in Section 6. 

2 Problem formulation 

2.1 Conventional formulation 

As mentioned, a mathematical formulation of the problem includes three groups of equations. Firstly 
we present them in the conventional form. 

fluid equations include the volume conservation law 

Ow 

— + divq + fl=0 (1) 

and the relation of the Poiseuille type obtained by integration of Navier-Stokes equations for a flow of 
viscous fluid in a narrow channel 

Q = —D(w,p)gia.dp. (2) 

Herein, w(x, t) is the channel width (fracture opening), q(x, t) is the flux vector through the fracture 
height, qi(x.,t) is the intensity of distributed sinks or sources (below this term will be assumed positive to 
account for leak-off), p(x, t) is the pressure, D is a function or operator, such that D(Q,p)gra,dp = 0, x 
denotes the vector of the position of a point on the surface of the flow, t is the time. The flux, divergence 
and gradient are defined in the tangent plane to the surface of the flow. 

Substitution of ([2]) into (p]) yields the lubrication (Reynolds) equation 

-^ - div (D(w,p) gradp) + qi = 0. (3) 

In hydraulic fracture problems, the opening is not known in advance. Thus its initial spatial distri- 
bution should be defined at start time io : 

w(x,t ) = wo{x), (4) 



where wo(x) is a prescribed function. 

The spatial operator in (|3]) is of the second order and elliptic. It requires only one boundary condition 
at the fluid contour L e . For instance, when neglecting the lag between the fluid front and the fracture 
contour, it may be the condition of the prescribed normal component q n of the flux: 

g„(x,t) = qo(x,t), x e L e , (5) 

where go( x , ^) is a known function at L e ; specifically, at the points of the fluid injection it is defined by 
the injection regime; at the points of the fluid front, coinciding with the fracture contour, we have w = 
and equation © implies go( x ,£) = 0. 

In conventional formulations (e.g. [H], [J, [27], [5], [5], [IS], 0, PU, QU, [ID]), to follow the fluid 
front propagation, authors use the equation of the global mass balance. Being a single equation, it do 
may serve for this purpose when considering 1-D problems with one point of the front to be traced. 
However in general, as emphasized in |19j - [22], it is preferable to employ the speed equation, which 
is formulated at each point of the fluid front. Even for 1-D problems, it provides advantages discussed 
in following sections. We shall not dwell on this issue here as the discussion below focuses on a 1-D 
problem. 

Solid mechanics equations define a dependence of the opening on the net pressure caused by defor- 
mation of rock: 

Aw = p, (6) 

with the condition of zero opening at each point x c of the fracture contour: 

w(x c ) = 0. (7) 

Commonly, the operator A in ([6]) is obtained by using the theory of linear static elasticity. As 
mentioned, when neglecting the lag, the condition of zero opening ([7]) replaces the condition of zero flux 
on the front. Henceforth, we shall consider this case and write x c = x* with the star marking that a 
quantity refers to the fluid front. 

Fracture mechanics equations define the critical state and the perspective direction of the fracture 
propagation. In the commonly considered case of the tensile mode of fracture, these are: 

Ki(x c ) = K IC , Kn(xc) = 0, (8) 

where Kj is the tensile stress intensity factor (SIF), Kjc is its critical value, Kjj is the shear SIF. 

The problem consists in solving the PDE ^ together with the elasticity equation ^ under the 
initial condition ((4]), boundary conditions (0, ([7]) and the fracture conditions ([8]). As mentioned, the 
global mass balance is usually employed instead of the speed equation to find a current position of the 
front. 

2.2 Modified formulation of fluid equations and boundary conditions 

The modified formulation |19| - |22| concerns mostly with the fluid equations and corresponding bound- 
ary conditions. It employs the primary quantity resulting from integration of the Navier-Stokes equations 
when a flow occurs in a narrow channel; this is the particle velocity. Its value v averaged across the 
channel height defines the flux q entering equations ([T]), ([2]), ([5]), because by definition 

q = lov. (9) 

Thus we may use the fluid particle velocity 

v=^ (10) 

w 



instead of the flux q. In terms of the particle velocity, the conservation law ([1} an d the Poiscuille type 
equation ([2]) become, respectively: 



dw 
~dt 



and 



+ div (ot) + qi = 0, (11) 



-— D(w,p) gradp. (12) 

w 



For the velocity component v n in a direction n, the equation (|12[) yields 

v n = --D(w,p)^. (13) 

w on 

In contrast with the flux, pressure and opening, the particle velocity is a smooth function near the 
fluid front. It follows from the fact that the particle velocity v(x,J equals to the front speed V* at each 
point x* of the front. In terms of the components «„», V*, normal to the front, we have: 

u„*(x*) = — ^ = V*(x*), (14) 

where x n * is the normal component of a point x* on the front, V* = |V*|. Hence, the particle velocity is 
finite at the front in common cases of the front propagation with a finite speed. Moreover, it is non-zero 
except for flows with stagnation points. 

The equation (|14p , where the normal component of the particle velocity is defined by (|13[) , presents 
the speed equation for the problem of hydraulic fracture: 

V*(x*) = --?—D(w,p)-^-. (15) 

Herein, n* is the unit normal to the front in the direction of its propagation at a point x*. Being the 
starting concept of the theory of propagating surfaces |28j , the speed equation is fundamental for proper 
tracing the hydraulic fracture propagation. 

The speed equation (|15p yields also important implications for numerical simulation of hydraulic 
fractures by finite differences (FD). Indeed, when at a time step we have known both x* and V*(x*), the 
equation (Ti"5|) becomes a boundary condition additional to the boundary condition (J7J) on the front. Thus, 
as noted in [19], a boundary value problem may appear overdetermined and ill-posed in the Hadamard 
sense |12j . To avoid difficulties, it is reasonable to use e-regularization, suggested and successfully used 

in rju, guj. 

The e-rcgularization is performed as follows. An exact boundary condition on the fluid front is 
changed to an approximate equality at a small distance r e behind the front. This approximate equality 
is obtained by combining the boundary condition at the fluid front, particular for a considered problem, 
with the speed equation, which is quite general. In practical calculations, the distance (absolute r e or 
relative e) is taken small enough to use the equality sign in the derived approximate condition. This 
gives us the e -regularized boundary condition near the front. The speed equation is also assumed to be 
met at the distance r E with an accepted accuracy. This gives us the e-regularized speed equation. The 
e-regularized boundary condition allows one to avoid the mentioned unfavorable computational effects; 
the e-regularized speed equation serves to find the front propagation. 

In this way, the boundary conditions (J7]) and (|TJ5) arc combined to obtain the e-regularized boundary 
condition [20], [21]: 

-D{w,p)dp = V*r e , (16) 

'PC, w 



where p 6 = p(r e ) is the pressure at the distance r e from the front. The e-rcgularized form of the speed 
equation (15) is: 

% {t) = ^ = -±D{ W ,p)^- . (17) 

at w on r E 

The equations (|16p and (jTTJ) actually employ the system moving with the front. Thus it is reasonable 
to re- write the lubrication equation (|11[) in this system. In it, the r-axis is directed opposite to the front 
velocity, while the other axis is tangent to the front. The connection between the temporal derivatives 
evaluated under constant x and r is given by the rule: 



— I = — I + K. — 

dt lx=corist dt \r=const * dr ' 



Then equation (fTTj) reads [20] : 



d In w dv n , , d In w 1 . „ 

-m—w + ^-^-ar-u*' (18) 

where using lnw serves to account for an arbitrary power asymptotic behavior of the opening 

w(r,t)=C(t)r a + 0(r 5 ), r -> 0, a > 0, 5 > a (19) 

near the front. The value of the exponent a is known in a number of important particular cases (see, 

e.g. [33], [5], PP, [H]), and <5 = 1 + a when the leak-off is neglected. 

For the asymptotics (19), it is reasonable, in addition to the particle velocity, to use the variable 

y(r, t) = [iu(x* — rn, t)] x / a , which is linear in r near the front. Finally the lubrication equation (18) near 

the fluid front becomes 

dy y dv n dy y 1 ' 

-TT7 = — tr + iVn-V*)-* 9l- (20) 

ot a or or a 

In 1-D cases, the equation (|20p is applicable to the entire fluid. In these cases, there is the only 

spatial coordinate x and it is reasonable to normalize x or, what is actually equivalent, r by the distance 

x*(£) from the inlet to the front. When using <^ = x/x*, the partial derivative evaluated under constant 

r is expressed via that under constant ? as: 

-I =-l +(1- )-- 

Then in terms of <; = x/x* = 1 — r/x*, the lubrication equation (|20|) in 1-D cases reads: 



dy = J_ 

dt x* 



Os a os 



y 1 -* - 



where we have omitted the subscript n in the notation of the particle velocity; tilda over a symbol 
marks that the corresponding function is considered to be a function of <;: y(s,t) — y(x*(l — s),t), 
v(<;,t) = w(a;*(l — ?),£), qi{s,t) = qi(x*{l — s),t). From now on, to simplify notation, we shall omit the 
tilda over functions depending on <;. Thus, the previous equation is written as: 



dy_ = l_ 

dt a;* 



dy y dv 

(?v; -v)- — 

ds a ds 



,i- 



■91, (21) 



o 



Note that when qi near the front decreases faster than w = y a , wc may divide (21) by y, obtaining 
the equation 

1 dy ?K - v dy 1 dv 1 

_ "a7 = "a a «*• (22) 

y dt x*y ds ai, o^ ay a 

In ([22]) . under the assumed asymptotics of qi, the term (dy/dt)/y, the factor (?V* — v)/(x*y) and the 
derivative dv/ds are finite at the fluid front. 



2.3 Nordgren problem in modified formulation 

The 1-D problem (Fig. 1) studied by Nordgren [25] is similar to that considered by Perkins and Kern 
[26] , improving the model of these authors by rigorous mathematical formulation, which includes finding 
the fracture length cc* as a part of the solution. Below we use the rigorous formulation by Nordgren and 
attribute it to this author. 

The fluid equations of the problem are (JTJ) - ((3j) with the operator D being the multiplier 



D(w,p) = hw 3 , 



(23) 



corresponding to the flow of Newtonian fluid in a narrow channel with an elliptic cross-section; for it 
ki = l/(7r 2 /z), where /i is the dynamic viscosity. 

The operator A in the solid mechanics equation (J5]) is also taken in the simplest form as the multiplier 
fc e = (2/wh)E/(l — v 2 ) in the linear dependence between the pressure and opening 



P 



(24) 



wcllborc 




elliptical crack 



Figure 1: Scheme of the Nordgren problem. 



Herein, h is the vertical length of a narrow elliptic channel (Fig. 1), E is the Young modulus, v is 
the Poisson's ratio of rock. For the dependence (|24|l , there is no need in a fracture criterion because it 
does not involve the fracture front. 

The initial condition Q in the 1-D case reads 



w(x,to) = w (x), 



(25) 



with wq(x) = ahead of the fluid front x*(to)- 

In view of ||3J), (|2"3"j) and (j2"4")h the boundary condition (|5|) of the prescribed influx go at the inlet x = 

becomes: 

kik e , , dw 3 1 , . . 



The boundary condition of zero opening involves the only point a;* = x*(t) of the fluid front: 

w(x*)=0. (27) 

The equation for the particle velocity (|12l) and the speed equation (TT5)) become, respectively, 

kik e dw 3 



3 dx 
and 



(28) 



_ dx^ _kik 1 dw 3 _, 
*~ dt ~ 3 dx l*=".W ^ 9J 

The exponent a in (19) equals 1/3 (see, e.g. [T7]) when the leak-off is neglected or even when it is 
described by the Carter's dependence and, consequently, singular at the crack tip. Then y = u> 3 , and 
in terms of the normalized independent variable £, the modified equation (21) for the 1-D problem (21) 
reads: 



dy = _}_ 

dt x* 



I \r \ d y q dv 



3y 2/3 qi, (30) 

where the particle velocity v is connected with y by equation following from 

k L k e dy 



3ic* 3? 



(31) 



In the new variables, the asymptotic equation (19) with a = 1/3 is written as y — C 3 (t)x*(l — q) + 
0((1 — ?) 1+K ) for <; — y 1, where n > for any leak-off tending to zero at the crack tip. Then (29) and 
(31) imply that near the front v = V* = ^pC 3 (£). Hence, the coefficient C 3 (t) is a multiple of the front 
speed: C 3 (i) = Trr-V*(t), and the asymptotics of y near the front is 

y(^ = ^-Mt)x,(t)(l-,) + 0((l-,) 1+K ), ?-►!. (32) 



In terms of the variables v, y and <;, the conditions (f25j) and (|26)) read, respectively: 

tf(?,«o)=lto(0, ( 33 ) 



W)K0,t) = »W, (34) 

where j/o(?) = Wg(<;a;*) for < ? < 1 and yo(<r) = ahead of the fluid front (? > 1). 

From (|32|) . ([29| it follows that the e- regularized boundary condition (fi~6|) and the e-regularized speed 

equation (|17|) arc, respectively: 

3 
y(t, 1 - e) = — — V*x*e, (35) 



kik e dy 
V*x* =- 



(36) 



3 a<r 

where e is a small relative distance from the front (e = r e /x*) and V* = dx*/dt. 

We need to solve the PDE (|30|) . where v is connected with y by equation (|3ip . under the initial 
condition (|33)l . the boundary condition (f34|) at the inlet and the e-regularizcd boundary condition (|35[) 
imposed at a small relative distance e from the fluid front. The regularized speed equation (|36j) serves 
to find the position of the fluid front. 

Emphasize that, considering BVP, we do not use the conditions ([2T]) and (|29p . not involving regular- 
ization, to avoid deterioration of the solution near the front. 



2.4 Self-similar formulation. Benchmark solutions 

As shown by Spcncc and Sharp [25], 1-D plane and axisymmetric problems may be reduced to a self- 
similar formulation in the case when there is no leak-off and the flux q at the inlet is proportional to 
a power ip(/3,t) = t" or exponential (p(f3,t) = e* 3 * function of time with constant f3. In the particular 
case of constant influx, j3 — 0. Representing a solution in the form w(t, <;) = (p(/3 w ,t)W(^), p(t, ?) = 
<p(/3 p ,t)P(<;) with separated temporal ip(j,t) and spatial <; = x/x* variables leads to equations with the 
only independent variable <; and with x* = Btp(P„,t). For an axisymmetric problem, x* is a current 
radius of the fracture. The constants W , /3_ and /3„ depend on a particular ID problem and /3. 

Actually, Nordgren employed this option for the case of constant influx (f3 — 0) ((22], Appendix 
C). We shall use the separation of variables with (p(/3,t) = t 13 to find benchmark solutions needed for 
further discussion. We include the case of non-zero leak-off by representing the leak-off term in separated 
variables, as well. 

In this way, for the considered problem, the speed equation yields /3„ = {3f3 w + 1)/2, B = yV(l)//3„, 
while the PDE (|30]) becomes the ordinary differential equation (ODE): 

where V(<;) is the self-similar particle velocity defined as V(?) = — ;j ^r- with K(c) = PF 3 (<;). For 
consistency, the leak-off term qi entering (|30)) is also taken as a function in separated variables qi = 

The condition (|34|) at the inlet defines the dependence of /^^ on the exponent (3 in the influx prescribed 
by q = Af 3 , where A is a constant: ^ = (1 + 2/3)/5. Then /3„ = (3/3 + 4)/5. Noting that the front 
propagates with the speed V*(t) = Bfi^*' 1 , this implies that the propagation speed is constant when 
/3„ = 1, what corresponds to j3 = 1/3. For a constant influx (/? = 0), we have the Nordgren's results: 
j3 w = 1/5, /3„ = 4/5 and the propagation speed changes proportionally to i~ 1//5 . 

2.4.1 Benchmark solution for zero leak-off 

For zero leak-off q\ = 0, (or Q/(?) = 0) and constant influx (j3 = 0, /3 M , = 1/5, /3„ = 4/5) the self-similar 
formulation serves to obtain the analytical solution |22j . In terms of physical quantities it is: 



«;(?, t) = w„ ^/Y(x/x^)t 1 / b , y = y n Y(x/x*)t 3 / 5 , p = k e w, 

(38) 

V = VnV^ix/x^t- 1 / 5 , V t (t) = 0.8^Vnt- 1/5 , X, (t) = ^3!„t 4 / 5 . 

Herein, x„ = (fc;fc e /4) 1/5 g T / t„ , iu n = q n t n jx n , y n = w^, v n = x n /t n , q n , and t n are normalizing 
length, opening, cubed opening, particle velocity, flux and time, respectively. The normalizing quantities 
q n , t n may be chosen as convenient. The dimcnsionlcss parameter £„ is defined by the prescribed influx 
<7o at the inlet [20] : 

£* = 1.3208446(q /gn) - 6 - 

Thus, when taking the influx go as the normalizing flux q n = <jo, one has £„ = 1.3208446. As above, 
? = x/x* is the relative distance from the inlet. The functions i></,(?) and Y"(<?) are given by the series: 

j—oo j—oo . 

MO = 0.^X) & i( 1 -^' nO = 0.6^]T^I(l-^, (39) 



where 6o = 1, &i = —1/16, and the next coefficients are found recurrently as 

b; 



'Jj+i 



3j + 4 



4j + l ^ 3j -_2fc + 6 , 

&j + 2^ 1 h-ibj-k+2 



40- + 1) 



1,2,. 



!) 



The series (39) rapidly converge. Five first terms with coefficients 60 = 1, 61 = —1/16, 62 = 
— (15/224)6i, 63 = —(3/80)62, ^4 = —(11/5824)63 provide the accuracy of seven significant digits in the 
entire interval of flow. The corresponding relative error is of order 10~ 5 even near the fluid front where 
F(?) — > 0. In further calculations, we use seven terms what guaranties that the relative error is of order 
10~ 7 . For £„ = 1, the normalized self-similar particle velocity v^,(<;) = V(?) and the cubed normalized 
opening Y(s) present the solution of the equation (37), when the variables in the latter are normalized 
in accordance with (38). Then Y(l) = 0, V(l) = 0.8, V(s) = -%% and Q;(?) = 0. Below we shall 
use the solution (38), (39) with q n = q a = 1, t n = 1, ki = 1, k e = 1; then x n = 4 -1 / 5 , w n = 4 1 / 5 , 
£* = 1.3208446. We shall call it the benchmark solution I. 

2.4.2 Benchmark solution for non-zero leak-off 

In this case, we prescribe the function W(?) and define the corresponding leak-off term Qi(<;) entering 
(j3"7| as such, for which the lubrication equation is satisfied by W(<;). Specifically, for a prescribed function 

W(?), the latter satisfies ([37]), when assuming Q;(<r) = -Y 1 / 3 [^ + VM ~^ y(1) ^ + ftj . The initial 
condition (j33j) becomes t/(?,to) = *o ™^ 3 W- We specify the function W(<r) by the expression: 

W r (?)=«(l-0 1/8 [1 + *W]. (4°) 

where w is a constant, s(?) = CV(1 — <r) + 0((1 — <r) 2 ) as <r — 5- 1, and C\y is a constant chosen in such 
a way that the particle velocity and leak-off term are positive in the entire flow region. The choice of 
W(<r) in the form (40) guaranties the asymptotic behavior (19) of the crack opening. 
The benchmark solution, corresponding to the choice (40), is: 

w{s,t) = W{<:)iP«, y = w\ p = k e w, B W = (1 + 2B)/5 1 0„ = (3/0 + 4)/5, 



, . 3k e kiui 3 (3B + 4) 3/3-1 /!,. , .. , . .„ . .. ,. ,\ 

«(?,*) = V 5 *^~ (3 [1 + S(?)1 " (1 " ?)[1 + S(?)]S W J ' 

go(t) = _ ^W + 4) ^ (1 + s(0)] 3 {s , (0) _ (1 + s(0)]/3}; 



k P kiuj 3 (3B + 4) 3/3-1 5k e kiu> 3 R 



23-4 

t^V- Z3/3 + 4 
91 ~ 3VK 2 (0 ^^5^" 



'aiF 3 dW 3 , \ dW 3 d 2 W 3 



\ a<r <9^ '«— iy ^ a<^ 



6/? + 3 



ir 



The initial condition (25) reads w(x, to) — W{x/x*)to , For certainty, in further calculations we set 
to = 1, k e = 1, ki = 1, w = 1 and consider the case of constant influx: j3 = (B w — 1/5, /3„ = 4/5). 
Below we shall use two choices of the function s(?) with the same Cw = — l/(96e). One of them 
s(?) = — (96e) (1 — q) corresponds to small leak-off; it will be referred as the benchmark solution II. 
The other s(<r) = — (96e) _1 (l — <;) +0.05(1 — <;) 2 describes notable leak-off; it will be referred as the 
benchmark solution III. 

In order to compare the benchmark solutions I, II and III, we evaluate the total fluid loss per unit 
time at the initial moment to = 1: Qi — L (&(?, l)cfe. For the benchmark solution I, we have Qi = 0; 
for the benchmark solution II, the total loss is Qi = 3.3 • 10~ 3 , Qi/qo = 4 • 10~ 3 ; for the benchmark 
solution III, the total loss becomes an order greater being Qi = 4.8 • 10~ 2 , Qi/qo = 4 • 10~ 2 . They also 
correspond to different distributions of the particle velocity in the flow region what has strong impact 
on the accuracy of calculations. We may characterize the variation of the velocity by the parameter 
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7„ = [max(u(<r,i)) -min(u(?,t)] J vz(<r,£)<*r 

Its values are: 7^ = 0.06 for the benchmark solution I, -f v = 0.02 for the benchmark solution II, 
and j v — 0.4 for the benchmark solution III. Thus, the velocity distribution is almost uniform for the 
benchmark solution II, and it is strongly non-uniform for the benchmark solution III; the benchmark 
solution I with zero leak-off presents an intermediate case. 

3 Spatial discretization. Stiffness analysis 

Henceforth, in accordance with e-regularization, we consider the spatial interval [0,1 — e] instead of 
[0,1] to avoid computational difficulties disclosed in [15] . [20] . Normally we shall set s = 10 -4 what 
guaranties that the relative error of the cubed opening is of order 10 -4 even near the fluid front. The 
spatial discretization of the problem employs representation of the interval [0, 1 — e] with N — 1 segments 
of the equal length h e = (1 - e)/(N - 1). The nodal points are ?j = j(l - e)/(N - 1) (j = 1, 2, ..., N). 
Then applying a finite difference approximation (say, the left-hand side approximation) to the spatial 
derivative(s) of the unknown function y, the PDE (21) yields a system of ODE in the vector Y of nodal 
values yj. 

3.1 Iterations in particle velocity 

Let us employ the fact that the particle velocity v, being the ratio of the flux and opening, which 
both decrease when approaching the fluid front, changes significantly less than these quantities. (As 
established in [22] . the particle velocity is almost constant in the entire flow region in problems of 
Nordgren and Spence & Sharp when neglecting the lag). Let us assume as a rough estimate v = V*. 
Then the ODE has the linear form: 



dY 
~~dt~ 



where 



A(t)Y + B(t), t>t , (41) 



A(t) =— DE, (42) 



D and E are diagonal and two-diagonal matrices, whose non-zero entries are defined, respectively, as: 
da = 1 — Ei, (i = 1, ..., N) and e^ = 1 for i = j, e^ = —1 for i = j + 1, (i, j = 1, ..., N). The eigenvalues 
Xj of the matrix A(i), defined by (42), are Xj = — djj "S.l , (j = 1,2..., N). All of them are negative, 



and the stiffness (condition) ratio Ua = "inf-AO 1 ' characterizing the stiffness of the system of ODE (41), 



Aj ui liiu iiiauix ±i.yi,j, uuimtjLi uy ^i^J, aie /\j — — Ujj m , 

max( — Xj ) 
min( — Xj) 

is 

rf ll(g) 1 M 

cLnn{£) £ 

Actually, it is possible to set e = when obtaining the system of ODE (41), because the initial (Cauchy) 
problem of solving it under the condition of zero opening Y(l) = is well posed. Still, effectively, 
when employing the condition Y(l) = 0, we arrive at the same system (41) with N — 1 unknowns and 
£ = 1/(N - 1). Then (43) becomes 

k a = N - 1. 

In a general case, the particle velocity is not constant along the fracture. Nevertheless, the stiffness 
ratio for large N can be estimated as 

K A ~-CA, (44) 
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where ca = ca (t) depends on the distribution of the velocity along the fracture and can be computed as 



dv , 



c A = v(0,t)\V*{t) + —(l,t) 

For the benchmark solutions under consideration, we have the stiffness ratio (|44[) independent of time t: 
c A = 0.52, c A = 0.96, c A = 1.43. Finally, employing the condition Y(l) = 0, the stiffness ratio is 
linear in N. 

Linear dependence on N is more favorable for solving ODE, than quadratic or cubic dependencies, 
which appear in other approaches. Even for N ~ 10 5 , the stiffness is of order 10 5 , what is quite acceptable 
in practical calculations. This indicates that solving the problem by iterations with the particle velocity, 
found at the stage of temporal integration, may be reasonable. However, to the moment, it is unclear 
how to properly organise the iterations to meet the BC of prescribed influx at the inlet. 

3.2 Spatial discretization with reduction to dynamic system of ODE 

Another way of solving the problem may consist in substitution of the velocity, defined by (31), into 
(30) and considering the PDF with the second partial derivative. The substitution may be employed 
cither for all terms including the particle velocity, or only for the term dv/d<;. Herein we employ the 
first option and obtain: 



dy_ _ kike_ ( d?y_ 1 
dt ~ xl \ V d<; 2 3 



dy dy 



<V| o 2/3, 



f q \-W /A m- (45) 



We compliment (45) with the regularized SE (36), which, as mentioned, serves to find the fracture 
length x*. When employing (45), it is extremely beneficial to re- write the SE (36) in the form, which 
similarly to (45) contains x\: 

dxl n kik e dy , .,„. 

^r = ~ 2 — a?U- e - (46) 

Then after a spatial discretization and using the boundary conditions at the inlet (34) and on the front 
(35), we arrive at a system of ODE in unknown values at nodes inside the flow region and the additional 
unknown x\. The initial conditions for the system are given by (33) and x\{t§) = xj . This opens 
the possibility to utilize well-developed methods for solving systems of ODE. Specifically, the Runge- 
Kutta method become available. Below we extensively use the new option for studying the accuracy 
and sensitivity. Emphasize that it has appeared only due to employing the (local) speed equation as 
the basis for tracing the front propagation. There is no such an opportunity when tracing the front in 
conventional ways (e.g. |16j . [2]) by using the global mass balance. Below we shall also see that the 
accuracy of calculations, based on the SE in the form (46), is significantly (more than an order) better 
than that obtained with the global mass balance. 

We may estimate the stiffness of the dynamic system (41), corresponding to (45) after spatial dis- 
cretization. The system is non-linear but taking into account the asymptotic behavior (32) and the 
leading term of the equation (45), we obtain the following expression for the matrix A(t): 

where the matrices D and E are defined above. In this case, it is possible to evaluate the condition 
number fcoo = ||A||oo ■ HA -1 ^ in the space l^ by using results from [13]. When taking e = 1/(N — 1), 
the estimation for large N is: 

fcoo = 2N 2 {lnN + ~f-l) + 0(N), N^oo, 
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where 7 is the Euler-Mascheroni constant. 

The dependence is close to quadratic in N what means much stronger growth of stiffness in the 
considered approach than in the previous one. Again, in a general case, there is no analytical formula 
for the stiffness ratio. It may be estimated as 



k a = c a N 2 + 0{N), N^oo, 



(47) 



where ca is a constant depending on a particular problem. We calculated the condition number k^ and 
the stiffness ratio ka numerically at to = 1 for the benchmark solutions I, II and III with e = 1/(N — I). 
The dependencies of koo and Ua on N for the benchmark solutions II and III are presented in Fig. [5J 
The graphs for the benchmark I are indistinguishable from those for the benchmark II. It can be seen 
that fcoo and kA comply with the asymptotic estimation (47). 
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Figure 2: Condition number Koo(N) and stiffness ratio ka{N) for the benchmark solution II (lines 
without markers) and III (lines with markers). 

Note that the quadratic dependence (47) on the number of nodes N follows from the proportionality 
of the pressure to the opening accepted in the Nordgren model. In the case when the pressure and the 
opening are connected by the exact equations of the elasticity theory, the dependence becomes cubic 
(see, e.g. [5]). For N = 10 5 , the stiffness ratio becomes of order 10 15 what is critical for most of the 
solvers. In this respect, employing iterations in the fluid velocity may be a reasonable strategy despite it 
requires repeated solving of ODE and accurate evaluation of the velocity after an iteration. Still, using 
the system (4-5), (46) looks beneficial in a vicinity of the fluid front in the general case of a 2D fracture. 
Then the number of nodal points may be taken small enough to avoid too stiff system. 

4 Accuracy of computations 

We shall use two mentioned approaches for approximation of the spatial derivatives. 

The first one, leading to the dynamic system (45), (46), has been suggested and explained when 
analyzing stiffness. It results in a well-posed initial (Cauchy) problem for the system of ODE what opens 
the possibility to solve the problem by well-established methods, like those of Runge-Kutta. Making use 
of this opportunity, we applied the standard MATLAB solver odel5s. It is based on the Runge-Kutta 
method and employs automatic choice of the time step. 

The other option consists in substitution of the velocity, defined by (31), into (30) only for the term 
dv/dq. Then the equation contains the propagation speed V* and the particle velocity v in the multiplier 
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by the first spatial derivative: 



at 



kik e d 2 y 1 



-ll 



dy 



P + ^-»»i- 3 » 2/3 



tu- 



rn 



We compliment (48) with the regularized speed equation (36), which after integration in time and 
accounting for the asymptotics (32) may be written as 



. fci.fcf, 



' lim y -^Hdt. 

to *-+! 1 - * 



(49) 



This form of equations is convenient for employing an implicit time stepping method, for example, 
the classical Crank-Nicolson scheme. The boundary conditions on a time step are obtained by using the 
discretized equations (34) (at the inlet) and (35) (at the e-regularized opening condition). The resulting 
non-linear algebraic system is solved by successive iterations. At each of the iterations, we consider a 
linear system by taking fixed values of y, v, V* and x* in the non-linear terms. At an iteration, for fixed 
coefficients in front of the spatial derivatives and leak-off term and for fixed coefficients in non-linear 
boundary conditions, the algebraic system is tri-diagonal. It is efficiently solved by the sweep-method. 
At the end of an iteration, we obtain new nodal values of y, which serve to evaluate new values of V* and 
v by using (31), and new value of a;* by using (49). The non-linear terms are iterated within a time step 
until the difference of values obtained on successive iterations becomes less than a prescribed tolerance. 
For a sufficiently small time step, the values obtained on the previous time step may serve as acceptable 
approximations, then one iteration is usually sufficient to meet the tolerance. 
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Figure 3: Relative error of the crack opening for the benchmark solution II: a) computations by Runge- 
Kutta method for the system (45), (46) after spatial discretization, b) computations by Crank-Nicolson 
method for the system (48), (49) after temporal and spatial discretization. 

We compare the accuracy of the solution, obtained by the standard Rungc-Kutta MATLAB solver 
odeA5s for the dynamic system, with that, obtained by the Crank-Nicolson method with iterations on a 
time step for the discretized system (48), (49) and discretized boundary conditions (34) and (35). In both 
cases, we used the same second order approximations for the first and second spatial derivatives. The 
number N of nodal points uniformly spaced on the interval [0, 1 — e] and the rcgularization parameter 
e, were also the same: N = 100 and e = 10~ 4 . Time stepping in the second approach was taken similar 
to that generated by the solver odel5s for solving the dynamic system of the first approach. The time 
interval [1, 100] was also the same for both approaches. Thus the conditions for the comparison were 
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practically the same. Fig. [3] illustrates the comparative accuracy of the two approaches. It presents the 
relative error Sw of the opening obtained by the first (Fig. |3^i) and the second (Fig. [3]b) approaches for 
the benchmark solution II. We see that the both approaches provide high accuracy: the relative error 
does not exceed 1.8 • 10 -6 . Still, the Crank- Nicolson scheme provides two-order less error (Sw < 10~ 8 ), 
except for the time close to the initial moment (to = 1) and nodes closest to the crack tip (<; = 1). 
(Surely, the error of the approach, using the Crank-Nicolson scheme, for small time and near the tip, 
may be decreased to the level 10 -8 , as well, by decreasing the first time steps and increasing the number 
of iterations within a time step). 

The Crank-Nicolson scheme served us also to compare the accuracy of results, obtained by using (i) 
the SE (49), and (ii) the global mass balance. In the normalized variables, the latter has the form: 

/*1 /*1 pt pt pi 

x*(t) / iu(?,i)d? = a;* (to) / w(^,t )d<; + / qo(t)dt- / x*(t) / qi(q,t)dqdt. 



Jo JO Jt J to Jo 

The calculations show that, when employing the global mass balance, the maximal relative error Sw 
of the opening becomes 1.4 • 10~ 3 . It is three-order greater than the error, obtained when employing the 
SE. This can be explained by an additional error induced by numerical integration over the interval [0,1]. 
This implies that using the local SE instead of the global mass balance is beneficial for the accuracy even 
in ID problems. 

Finally we analyze the dependence of the accuracy on the distribution of the fluid velocity along the 
fracture. Table 1 presents the maximal relative error of the opening Sw and the fracture length fa* for 
the benchmark solutions I, II and III. The data are obtained by using the first approach. The second 
line of the table contains the maximal relative variation "f v of the fluid velocity. It can be seen that the 
variation notably influences the accuracy. Even in the cases of the benchmark solutions I and II, when the 
variation j v is of the same order, there is significant (an order) difference in the accuracy. For strongly 
non-uniform distribution, corresponding to the benchmark solution III, the relative error is three orders 
greater than that for the most uniform distribution, corresponding to the benchmark solution II. This 
implies that numerical modeling of fractures with notable variations of the fluid velocity requires tests 
with growing density of the spatial mesh. 





benchmark I 


benchmark II 


benchmark III 


lv 


0.06 


0.02 


0.4 


Sw 


4.86 -10- 5 


1.54 -10 -6 


5.04- 10~ 3 


fa* 


7.23 -10- 5 


1.97 -10- 6 


6.84 • 10- 3 



Table 1: Accuracy of the crack opening Sw and crack length fa* for various benchmark solutions. 

5 Sensitivity analysis 

For zero leak-off, the only parameter, defining the solution, is the influx at the inlet qg. Consider its 
perturbed value go + Ago, so that the relative change of the influx is Sq = Aq n /q . We are interested 
in finding relative changes of the opening Sw = Aw/w, the pressure 5p = Ap/p, the fracture length 
fa* = A.t*/.t*, the particle velocity Sv = Av/v and the front speed SV = A14/V*. From the analytical 
solution (38), (39), we easily obtain: 

Sw = Sp = (1 + 5£*) 2/3 - 1, Sv = 8V* = fa* = S^. 
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In Sec. 2, it was stated that £„ is proportional to q . Hence, 1 + <5£„ = (1 + <5go) 3 ^ 5 , and finally the 
relative changes are: 

Sw = Sp=(l + Sq ) 2/5 - 1, Sv = SV„ = 5x* = (1 + Sq Q ) 3/5 - 1. 

For a small relative change of the influx, to the accuracy of terms of order O (Sq^) , we have: 

5w = Sp = 0A5q o - 0.126qo, 5v = 8V* = Sx* = 0.6Sq - 0.12^. (50) 

This implies that the relative changes equal approximately to one-half of the relative change of the 
influx. Note that the disturbances of the solution change its sign when 5qo changes the sign. 

In cases, when the disturbance 5qo oscillates in time with the amplitude A and the angular frequency 
(jj q as Sqo = Asm(ui q (t — to)), we may numerically evaluate its influence on the solution by employing 
any of the discussed approaches. The results, obtained when solving the discretized dynamic system 
(45), (46) by using standard MATLAB solver odel5s are presented in Fig. |4j 
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Figure 4: Relative fluctuations of the crack opening Swq and the crack length 5x*, caused by the periodic 
disturbance of the influx: a) benchmark solution II, b) benchmark solution III. 

The graphs are plotted for the relative amplitude of perturbation A/qo = ±0.1 and uj q = n/4. They 
present the relative fluctuations of the crack opening, Swo = Aw(0,t)/w(Q,t), at the inlet point (<j = 0) 
and relative fracture length 5x* in time for the benchmark solution II (Fig. |4^i) and III (Fig. [4J5) . Solid 
lines refer to the case when A/qo = +0.1, dashed lines to A/qo = —0.1. It can be seen that for the 
periodic perturbations the crack opening Swo is an order less than the amplitude A of the perturbation 
Sqo itself, while the position of the crack, fe*, is two orders less. It is worth noting that the fluctuations 
depend on the sign of the amplitude A, and the change of the fracture length tends to a constant 
value with growing time. It can be seen that, when the amplitude A changes its sign, the limit of the 
fluctuation of the fracture length for small leak-off (benchmark II) also changes its sign, while for the 
large leak-off (benchmark III) the sign is the same. 

6 Conclusions 

The results presented show that the modified formulation of the problem [19] -- [22], employing the 
suggested variables, speed equation and e-regularization, extends the opportunities for better numerical 
modeling of the hydraulic fracture propagation. Specifically, the improvements may include: 
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(i) drastic decrease of stiffness of ODE, obtained after spatial discretization, by employing iterations 
in the particle velocity; 

(ii) avoiding iterations in non-linear terms and using highly efficient standard solvers by including 
the SE into the system of ODE as suggested in Sec. 3; 

(iii) increasing the accuracy of time-stepping schemes of Crank-Nicolson type by using the SE instead 
of the commonly used global mass balance. 

The advantages of the modified formulation have been illustrated by considering the simplest ID 
problem by Nordgren. Meanwhile, the suggested formulation and approaches, being general, they are 
applicable to 2D fractures. They are of special interest for modeling the area behind the fluid front, 
where gradients of the pressure and opening arc high. 
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